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Abstract 

Spike-Timing Dependent Plasticity (STDP) is characterized by a wide range of temporal kernels. However, much of the 
theoretical work has focused on a specific kernel - the "temporally asymmetric Hebbian" learning rules. Previous studies 
linked excitatory STDP to positive feedback that can account for the emergence of response selectivity. Inhibitory plasticity 
was associated with negative feedback that can balance the excitatory and inhibitory inputs. Here we study the possible 
computational role of the temporal structure of the STDP. We represent the STDP as a superposition of two processes: 
potentiation and depression. This allows us to model a wide range of experimentally observed STDP kernels, from Hebbian 
to anti-Hebbian, by varying a single parameter. We investigate STDP dynamics of a single excitatory or inhibitory synapse in 
purely feed-forward architecture. We derive a mean-field-Fokker-Planck dynamics for the synaptic weight and analyze the 
effect of STDP structure on the fixed points of the mean field dynamics. We find a phase transition along the Hebbian to 
anti-Hebbian parameter from a phase that is characterized by a unimodal distribution of the synaptic weight, in which the 
STDP dynamics is governed by negative feedback, to a phase with positive feedback characterized by a bimodal 
distribution. The critical point of this transition depends on general properties of the STDP dynamics and not on the fine 
details. Namely, the dynamics is affected by the pre-post correlations only via a single number that quantifies its overlap 
with the STDP kernel. We find that by manipulating the STDP temporal kernel, negative feedback can be induced in 
excitatory synapses and positive feedback in inhibitory. Moreover, there is an exact symmetry between inhibitory and 
excitatory plasticity, i.e., for every STDP rule of inhibitory synapse there exists an STDP rule for excitatory synapse, such that 
their dynamics is identical. 
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Introduction 

Spike timing dependent plasticity (STDP) is a generalization of 
the celebrated Hebb postulate that "neurons that fire together wire 
together" to the temporal domain, according to the temporal 
order of the presynaptic and postsynaptic spike times. A 
temporally asymmetric Hebbian (TAH) plasticity rule has been 
reported in experimental STDP studies of excitatory synapses [1— 
3], in which an excitatory synapse undergoes long-term potenti- 
ation when presynaptic firing precedes the postsynaptic firing and 
long-term depression is induced when the temporal firing order is 
reversed, e.g., Figure 1A. 

Many theoretical studies [4-9] that followed these experiments 
used an exponentially decaying function to represent the temporal 
structure of the STDP. Throughout this paper we term this STDP 
pattern the "standard exponential TAH". Giitig and colleagues 
[7] also provided a convenient mathematical description for the 
dependence of STDP on the synaptic weight in the standard 
exponential TAH STDP rule: 



Aw=+Xf+(w)K(At) (1) 
is:(A0=e HA ' l/T (2) 



f + (w) = (l-w)" (3) 



/_(«■) = aw" (4) 

where we[0,l] is the dynamic parameter that describes the 
synaptic strength; Aw is the modification of w following pre (— ) or 
post (+) synaptic firing; At is the time difference between the 
presynaptic and postsynaptic firing; X is the learning rate; T is the 
temporal decay constant and //e[0,l] and a>0 are dimensionless 
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Figure 1. Illustration of different STDP temporal kernels (A' , ) 
as defined by equations (7) and (8) with the "standard 
exponential TAH" as a reference. Each plot (normalized to a 
maximal value of 1 in the LTP branch) qualitatively corresponds to some 
experimental data. In all plots, the blue curve represents the 
potentiation branch K + , the red curve represents the depression 
branch — K_ and the dashed black curve represents the superposition/ 
sum of K + —K-. For simplicity, all plots were drawn with the same 
x = 20ms. (A) The "standard exponential TAH" [1,18]. (B) ±K ± (At; 1,t) 
Alternate approximation to the standard exponential TAH [1,18], (C) 
±K ± (At; — 1,t) Temporally asymmetric Anti-Hebbian STDP [15]. (D) 
±K ± (Ar;0.15,r) TAH variation [12,19]. (E) +K ± (A?;0,t) Temporally 
symmetric Hebbian STDP [16,17]. (F) ±K ± (At; -0.75,t) Variation to a 
temporally asymmetric Anti-Hebbian STDP [19] 
doi:10.1371/journal.pone.0101109.g001 

parameters of the model that characterize the weight dependent 
component of the STDP rule. This representation introduces a 
convenient separation of variables, in which the synaptic update is 
given as a product of two functions. One function is the temporal 
kernel of the STDP rule, i.e. K(At), and the other is the weight 
dependent STDP component, i.e. f+(w). For convenience, 
throughout this paper we shall adopt the notation of Giitig and 
colleagues for the weight dependence of the STDP rule,/ + /_ (ir), 
equations (3) - (4). This function, / + /_(vf), is characterized by two 
parameters: the relative strength of depression - a, and the degree 
of non-linearity in w of the learning rule - fl. Note, that other 
choices for / + /_ (w) have also been used in the past [5],[10],[1 1]. 

Properties of the "standard exponential TAH" 

As previously shown [6,7], the standard exponential TAH 
model can generate positive feedback that induces bi-stability in 
the learning dynamics of an excitatory synapse. For a qualitative 
intuition into this phenomenon, consider the case of a weight- 
independent STDP rule, also termed the additive model, i.e., 
[1 = 0. If the synaptic weight is sufficiently strong, there is a 
relatively high probability that a presynaptic spike will be followed 



by a postsynaptic spike. Hence, causal events (i.e., A/>0 post 
firing after pre) are more likely to occur than a-causal events (with 
At<0). Because the STDP rule of the standard exponential TAH 
model implies LTP for A?>0 there is a greater likelihood for LTP 
than for LTD. Thus, a "strong" synapse will tend to become 
stronger. On the other hand, if the synaptic weight is sufficiendy 
weak, then pre and post firing will be approximately uncorrelated. 
As a result, the stochastic learning dynamics will randomly sample 
the area under the STDP temporal kernel. Here we need to 
consider two types of parameter settings. If the area under the 
causal branch in equation (1) is larger than the area under the a- 
causal branch, a < 1 , the net effect is LTP for weak synapses as 
well. Thus, in this case, all synapses will potentiate until they reach 
their upper saturation bound at 1 . Hence, the regime of a < 1 , in 
this case, is not interesting. On the other hand, if the area under 
the a-causal branch is larger than the area under the causal 
branch, cc > 1 , random sampling of the STDP temporal kernel by 
the stochastic learning dynamics (in the limit of weakly correlated 
pre-post firing, mentioned above) will result in LTD. Thus, in the 
interesting regime, (X>1, a "weak" synapse will tend to become 
weaker; thus producing the positive feedback mechanism that can 
generate bi-stability. 

It was further shown [7] that this positive feedback can be 
weakened by introducing the weight dependent STDP component 
via the non-linearity parameter fi in equations (3) and (4). Setting 
H>0 decreases the potentiation close to the upper saturation 
bound and decreases the depression close to the lower saturation 
bound; thus, for sufficiendy large values of f-i the learning dynamics 
will lose its bi-stability. 

Experimental studies have found that the temporally asymmet- 
ric Hebbian rule is not limited to excitatory synapses and has been 
reported in inhibitory synapses as well [12]. Similar reasoning 
shows that in the case of inhibitory synapses the standard 
exponential TAH induces negative feedback to the STDP 
dynamics. It was shown [13] that this negative feedback acts as 
a homeostatic mechanism that can balance feed-forward inhibi- 
tory and excitatory inputs. Interestingly, Vogels and colleagues 
[14] studied a temporally symmetric STDP rule for inhibitory 
synapses, and reported that this type of plasticity rule also results in 
negative feedback that can balance the feed-forward excitation. 
This raises the question whether inhibitory plasticity always results 
in a negative feedback regardless of the temporal structure of the 
STDP rule? On the other hand, theoretical studies have shown 
that the inherent positive feedback of excitatory STDP causes the 
learned excitatory weights to be sensitive to the correlation 
structure of the pre-synaptic excitatory population - for different 
choices of STDP rules [5,7,10]. Does STDP dynamics of 
excitatory synapses always characterized by a positive feedback? 

Outline 

Although theoretical research has emphasized the standard 
exponential model, empirical findings report a wide range of 
temporal kernels for both excitatory and inhibitory STDP; e.g., 
[1,12,15-19], (see also the comprehensive reviews by Caporale 
and Dan [20] and Vogels and colleagues [21]). Here we study the 
effect of the temporal structure of the STDP kernel on the 
resultant synaptic weight for both excitatory and inhibitory 
synapses. This is done in the framework of learning of a single 
synapse in a purely feed-forward architecture, as depicted in 
Figure 2. First, we suggest a useful STDP model that qualitatively 
captures these diverse empirical findings. Below we define our 
STDP model. This model serves to study a large family of STDP 
learning rules. We derive a mean field Fokker-Planck approxima- 
tion to the learning dynamics and show that it is governed by two 
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critical for the analysis below. Other types of functions may serve 
as well. The "Skew-Normal distribution" is defined by: 



SN(At; £,t, 9 



1 



x\>2n 



-e 2 



m 2 



l+erff^ 



(6) 



where £ is the temporal shift, x is the temporal decay constant, and 
<f is a dimensionless constant that affects the skewness of the curve 
and erf (x) is the Gaussian error function. It is also useful to reduce 
the number of parameters that define the STDP temporal kernel. 
Thus, we define: 



K + {At;d,x) = SN(At;2T8(\- 



9 2 ),T,20fJ 3 ) 



(7) 



Figure 2. Model architecture. The STDP dynamics of a single either 
excitatory or inhibitory synapse is studied in purely feed-forward model. 
In all of the simulations presented here, the activity of the presynaptic 
inputs is modeled by a homogeneous Poisson process, with mean firing 
rate r pre = lOspikes/s. The synaptic weights of all synapses except one is 
kept fixed at a value of 0.5. The post synaptic neuron is simulated using 
an integrate and fire model as elaborated. See Methods for further 
details. 

doi:10.1371/journal.pone.0101109.g002 



global constants that characterize the STDP temporal kernel. 
Stability analysis of the mean-field solution reveals that the STDP 
temporal kernels can be classified into two distinct types: Class-I, 
which is always mono-stable, and Class-II that can bifurcate to bi- 
stability. Finally, we discuss the symmetry between inhibitory and 
excitatory STDP dynamics. 

Results 

Generalization of the STDP rule 

In order to analyze various families of STDP temporal kernels 
found in experimental studies [1,12,15-19] we represent the 
STDP as the sum of two independent processes: one for 
potentiation and the other for depression. The synaptic update 
rule that we use throughout this paper is given by: 



Aw = X\f+ (w)K + (At) -/_ (w)K_ (At) 



(5) 



Note that the main distinction between equation (5) and (1), is 
that here, equation (5), the + / — signs denote potentiation and 
depression, respectively; while in equation (1) the +/— signs 
denote causal/a-causal branch. Thus, in our model for every At 
the synapse is affected by both potentiation and depression; 
whereas, according to the model of equation (1) the synapse 
undergoes either potentiation or depression - depending on the 
sign of At. For example, in the standard exponential TAH model 
defined above in equation (2), this generalization implies a 
K ± (At) = ®(±At)e + A 'l T temporal kernel, where &(At) is the 
Heaviside step function. For the weight dependent STDP 
component (f+(w)) we adopt the formalism of equations (3) and 
(4). Below we describe a workable parameterization for the STDP 
temporal kernel. 

Skew-Normal kernel 

Here we used the Skew-Normal distribution function to fit the 
temporal kernels of the STDP rule, K+ (At). Note that the specific 
choice of the Skew-Normal distribution is arbitrary and is not 



K_ (At; 6,%) =K + [At; —0,i \ 1 + 0.5(1 — 0 



(8) 



where 9e[ — 1,1] is a single continuous dimensionless parameter of 
the model that characterizes the STDP temporal kernel and i is 
the time constant of the exponential decay of the potentiation 
branch. The mapping of £cc0(l — 0 2 ) ensures that the temporal 
shift parameter, will be zero for 9 = — 1 ,0, 1 . In order to 
obtain temporally symmetric Mexican hat STDP rule for 



= 0 one needs to demand Tdep>' r ; 



L pot; 



where t 



and TH e p — T (l +0.5(1 -9 2 ) ). We also required T d , 



pot " 

^pot 



for 



9 = + 1, in order to be compatible with several previous studies. 
This reduction in parameters was chosen in order to capture the 
qualitative characteristics of various experimental data; however, 
other choices are also possible. Figure 1B-F illustrates how one can 
shift continuously from a temporally asymmetric Hebbian kernel 
(6=1 ,Figure 1 B) to a temporally asymmetric anti-Hebbian kernel 
(9 = — 1, Figure IF). Figure 1A shows the temporal kernel of the 
standard exponential TAH model, compare with 9=\, Figure IB. 

"Mean field" Fokker-Planck approximation 

We study the STDP dynamics of a single feed-forward synapse 
to a postsynaptic cell receiving other feed-forward inputs through 
synapses that are not plastic. We assume that all inputs to the cell 
obey Poisson process statistics with constant mean firing rate, r pre ; 
that the presynaptic firing of the studied synapse is uncorrelated 
with all other inputs to the postsynaptic neuron; and that the 
synaptic coupling of a single synapse is sufficiently weak. The 
STDP dynamics is governed by two factors: the STDP rule and 
the pre-post correlations. To define the dynamics one needs to 
describe how the pre-post correlations depend on the dynamical 
variable, it'. Under the above conditions one may assume that the 
contribution of a single pre-synaptic neuron that is uncorrelated 
with the rest of the feed-forward input to the post-synaptic neuron 
will be small. Thus, it is reasonable to approximate the pre-post 
correlation function (see Methods - equation (26)) up to a first 
order in the synaptic strength w (e.g., [8,22-24]), yielding: 



r(A) = <p.„(0p» ort (t+A)> = 



"pre* post 

(1+h 7 (A)),A>0 



' pre ' post? 



A<0 



(9) 



where p pre j post is the instantaneous firing of the pre/post synaptic 
cell represented by a train of delta functions at the neuron's spike 
times (see Methods), r pre jr posl is the pre /post synaptic mean firing 
rate; and the function y(A) describes the change in the conditional 
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mean firing rate of the postsynaptic neuron at time A + 1 following 
a presynaptic spike at time t. Note that we use upper case T to 
represent the full pre-post correlations, r = (.p pre P pos t) , whereas y 
denotes the first order term in the synaptic weight, w, of these 
correlations. 

In the limit of a slow learning rate, X— >0, one obtains the mean- 
field Fokker-Planck approximation to the stochastic STDP 
dynamic (see Methods - equation (27)), and using the linear 
approximation of the pre-post correlations, equation (9), yields: 



<vv>(w) =Xr pre r pos \f + (w) 



K+(A)dA + w 



y(A)K + (A)dA) 



(10) 



-/_(w)Q A"_(A)</A + wJ y(A)K_(A)dAj 
= ht pre r post [/+ (w) (K+ + wyK+ ) -/_ (w) (K- + wyK- ) ] 



whereF= J™^ F(A)dA denotes the mean over time (using equa- 
tion (9) with y(A) = 0 for A < 0). 

In our choice of parameterization, K + i _ are set to have the 
same integral; i.e., K + =K_ =K. The difference between the 
strength of potentiation and depression of the STDP rule is 
controlled by the parameter a. (equation (4)). Substituting 
expressions (3) & (4) into equation (10) yields: 

(wXw)=Xr pre r post K[(l-w)>'(l + wX + )-aw>'(l + wX^)} (11) 



where X + i _ =yK + i _ / K + / _ are constants that govern the 
mean-field dynamics. A fixed point solution, w*, of the mean- 
field Fokker-Planck dynamics, <tv>(w) =0, satisfies: 



l+w*X+ (\ 

Y+w*xl 



IV* 



(12) 



Numerical simulations - the steady state of STDP 
learning 

We performed a series of numerical simulations to test the 
approximation of the analytical result of the mean-field approx- 
imation at the limit of vanishing learning rate, using a conductance 
based integrate-and-fire postsynaptic neuron with Poisson feed- 
forward inputs (see Methods for details; a complete software 
package generating all the numerical results in this manuscript can 
be downloaded as File SI). We simulate a single postsynaptic 
neuron receiving feed-forward input from a population of 
A f g=120 excitatory neurons and Nj = 40 inhibitory neurons 
firing independendy according to a homogeneous Poisson process 
with rate r pre =lOspikes/s. All synapses except one (either 
excitatory or inhibitory) were set at a constant strength (of 0.5). 
The initial conditions for the plastic synapse were as specified 
bellow. 

We first estimated the spike triggered average (STA) firing rate 
of a single presynaptic neuron triggered on postsynaptic firing, in 
order to approximate the function y(A), equation (9). Figure 3 
shows the STAs of excitatory (A) and inhibitory (B) synapses for 
varying levels of synaptic weights (color coded), as were estimated 
numerically (dots). The dashed lines show smooth curve fits to the 
STA. The specific temporal structure of these curves depends upon 



particular details of the neuronal model. Nevertheless, the linear 
dependence on the synaptic weight is generic for weak synapses; 
thus, in line with the assumed linearity of the model, equation (9). 
The STA shows the conditional mean firing rate of the pre- 
synaptic neuron, given that the post fired at time t = 0. In the limit 
of weak coupling, vv— >0, pre and post firing are statistically 
independent and the conditional mean equals the mean firing rate 
of the pre, r pre = \Qspikes / s. For an excitatory synapse, as the 
synaptic weight is increased the probability of a post spike 
following pre will also increase. Consequently, so will the 
likelihood of finding a pre spike during a certain time interval 
preceding a post spike. Hence, the STA of an excitatory synapse is 
expected to show higher amplitude for stronger synapse (as shown 
in A). Correspondingly, the STA of an inhibitory synapse is 
expected to show a more negative amplitude for stronger synapse 
(as shown in B). 

To fit the STA with an analytic function we used y(A) = 

a y sin((a)yA^ y ^Je~ A /' c '>, with the fitted parameters a y ,(O y ,ij/ y ,Xy for 

both the inhibitory and excitatory cases. This ad hoc approxima- 
tion serves to enable the numerical integration that calculates the 
constants X + i_ that govern equation (11) for the mean field 
approximation to the STDP dynamics. 

All the richness of physiological details that characterize the 
response of the post-synaptic neuron affect the STDP dynamics 
only via the two constants X + and X-. These two constants 
X + /_ denote the overlap between the temporal structure of the 
pre-post correlations, y(A), and the temporal kernel, K + i_, of the 
potentiation/depression kernel, respectively, Figure 4. Conse- 
quently, as y(A), is positive for excitatory synapses and negative for 
inhibitory synapses - so are the constants X + /_. In addition, as 
the correlations in our model are causal, y(A)=0for t<0, the 
constant X + (X-) is expected to decay to zero when the STDP 
kernel K + (KJ) vanishes from the causal branch, 6-* 1 (6-* — 1). 
For the specific choice of parameters in our simulations, \X + \ 
obtains its maximal value at 0=1. However, one may imagine 
other choice of parameters in which \X + | will obtain its maximal 
value at 0e(O,l). Note, from Figure 4, that the crossing of the X + 
and X_ curves, is coincidentally almost the same for both synapse 
types, and is obtained at « —0.2. The significance of this point is 
discussed below. 

Fixed point solutions for the STDP dynamics 

Figure 5 shows w* as a function of a. for different values of fi 
(color coded, note that a and ji are the parameters that 
characterize the synaptic weight dependence of the STDP rule, 
equations (3) and (4)). The panels depict different STDP setups 
that differ in terms of the temporal kernels as well as the type of 
synapse (excitatory/inhibitory). These two factors affect the mean 
field equations via X + i _ . The dashed lines show the solution to 
the fixed point equation, equation (12), using the numerically 
calculated X + /_. The fixed points were also estimated numeri- 
cally by directly simulating the STDP dynamics in a conductance 
based integrate and fire neuron (circles and error bars). 

For the estimation of the steady state value of the synaptic 
weight, the simulations were set to run for 5 hours of simulation 
time, which, according to manual offline analyses of convergence 
time scales, is much more than twice the time required for the 
system to converge and fluctuate around its steady state. The 
circles and error bars depict the mean ± standard deviation of 
the synaptic weight, as estimated from the last 2.5 hours of the 
simulation (weights were recorded at a 1 Hz sample rate). Note the 
high agreement between the fixed point solution (vv*) of equation 
(12), and the asymptotic synaptic weight (wq) as estimated by the 
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Figure 3. Spike Triggered Average (ST A) of a single presynap- 
tic input. The conditional mean firing rate of the presynaptic cell given 
that the postsynaptic cell has fired at time t — Q, is plotted as function of 
time. (A) Excitatory synapse (B) Inhibitory synapse. Each set of dots 
(color coded) is the conditional mean firing rate calculated over 
1000 hours of simulation time with fixed synaptic weights and 
presynaptic firing rates on all inputs. The different sets correspond to 
a different presynaptic weight [w*) on a single synapse on which the 
STA was measured. The respective dashed lines show the numerical 
fitting of the form STA(w x ,A) = r pre {\ + w*y(A)) where y(A) takes the 



revised formula: a,, sinl 



)e- A '-». 



u r e l a((a>yAy y \e A / T '. For every type of synapse, 

excitatory (in A) and inhibitory (in B), the parameters describing y(A), 
namely a y ,(Oy,\ji„,x y , were chosen to minimize the least square 
difference between the analytic expression and the numerical 
estimation of the STA. These parameters were then used to calculate 

X + l- 

doi:10.1371/journal.pone.0101109.g003 



The panels of Figure 5D show the results of the symmetric 
STDP with 0 = 0- note, that the dynamics of inhibitory synapse 
under the symmetric STDP rule, is characterized by a one to one 
function w* of a. corresponding to negative feedback, as previously 
reported [14]. 

Stability of the fixed point solution 

The stability of the fixed point solution w* to equation (12) is 
determined by the sign of the partial derivative of the 
dynamical equation, equation (11), with respect to the synaptic 
weight: 



a<w> 

+ ( 1 - w*fX+ - ot/iW> -1(1 + w*X_ ) - aw*nX_ ] 



(h>*) = Xr pre r post K [ - fi{ 1 - w*f- 1 ( 1 + w*X + ) 



On the other hand, examination of Figure 5 suggests that the 
stability of the fixed point is governed by the sign of dw* /dot. 
Taking the logarithm and the derivative with respect to a. of both 
sides of equation (12), one obtains: 



(14) 



8w* 


1 




x+ 


X_ fi ji 




8a 


a 


1 + 


w*X+ 1 


+w*X_ 1— w* w* 


sign 


~8w*~ 


= sign 


x+ 


X_ n 


8&L 


}+w*X+ 


l+w*X- l-w* 



(15) 



where the last equality holds as 0£>0. At the fixed point, 
substituting equation (12) into equation (13) one obtains: 



numerical simulation (regression coefficient of l + 5-10~ 4 with 
R 2 > 0.999 when performing a regression test on the entire set 
{wcH'*} presented in each of the panels). 

The panels of Figures 5A and 5B compare the standard 
exponential TAH rule of equation (2), in A, and our current 
STDP model with 6=1 in B, for a representative set of 
parameters {(cc,,^-)} applied to the examined synapse (middle 
column for inhibitory synapse and right column for excitatory 
synapse). Note that some lines may overlap each other near the 
boundaries: H'o = 0,l. As is evident from the figures, the results 
of the two models coincide. In particular, the Hebbian STDP 
dynamics of inhibitory synapses is characterized by a one to 
one function w* of cc and there is no bi-stability, as previously 
reported [13,14]. On the other hand, the Hebbian STDP of 
excitatory synapses is characterized by bi-stable solutions at 
low levels of ji below a certain critical value, see e.g., [6,7]. 
Thus, the current model with f?=l coincides with previous 
results. 

The panels of Figure 5F show the results of a temporally 
asymmetric /Irefe'-Hebbian STDP with 9 = — 1 . In striking contrast 
to the Hebbian STDP, in this case, inhibitory plasticity is 
characterized by bi-stability whereas, the excitatory plasticity is 
characterized by mono-stability. 

The panels of Figures 5C and 5E explore two other types of 
asymmetric rules (Hebbian and ^taz'-Hebbian respectively). These 
results show similar behavior as 5B and 5F in terms of the 
classification of STDP kernels discussed in the next section. 




-0.2 1 1 1 1 1 

-1 -0.5 0 0.5 1 

6 

Figure 4. Mean field constants A of equation (11) for the 
excitatory and inhibitory synapses of the neuronal model used 
in our numerical simulations, as a function of 0. These values 
were calculated using numerical integration (see File S1) with 
K + i_(At; 0,t) as defined by equations (7) and (8), with z = 20msec as 
set throughout the simulations, and with the fitted formula for y(A). 
doi:10.1371/joumal.pone.0101109.g004 



PLOS ONE | www.plosone.org 



5 



July 2014 | Volume 9 | Issue 7 | e101 109 



The Effect of STDP Temporal Kernel 



STDP kernel 



Inhibitory synapse 



Excitatory synapse 



AO 

st. TAH 



BO 

e = 1.00 



^7 



CO 

e = 0.75 




A 


e = - 




V 



E0 



F0 

e = -ioo 



3° 0.5 





\ 




~x„ 1 
--a 


s 


i s ~fc-.,i--i_-_-*--.,..S 




1 ^— «-—»—* 


■ 






- — — ^« — 


U 

i 

i 

1 ^4^-- Z. 




H = 0.00 ^V'i* 

n = o.oi ^Itf 

• n=o.o2 /^;-»». Ht 

H = 0.05 ,* V -1-^ 

n=o.io //L V_ 

|j = 0.2Q / .-i^ T *}-*-- T '* 



A1 i - . - ■ . 



C1 



F1 



A2 



//I I \I "1 

St ! ^ 



B1 



62 



.4' 




C2 



07 f-v*. 



D2 



XT'! 



£2 



A 



F2 



— 



0.7 0.8 0.9 



1 

a 



1.1 1.2 



0.7 0.8 0.9 1 1.1 1.2 

a 



Figure 5. The fixed point solution (ir ) of equation (12) (dotted lines), is compared to the asymptotic synaptic weight [wo) (circles), 
of a single synapse learning dynamics for various learning rules as defined by equation (5). Each of the panels in the middle column (for 
inhibitory synapse) and in the right column (for excitatory synapse) explores the weight dependent STDP component, /+ (w) of equations (3) and (4), 
for representative set of /( (shown by different colors as depicted in the legend) as a function of ct. The different rows correspond to different STDP 
kernels, K±(At) as shown by the panels in the left column. The circles and error bars represent the mean and standard deviation of the synaptic 
weight [wq), calculated over the trailing 50% of each learning dynamics simulation (see Methods). The mean field constants {X + ,X-} were 
numerically calculated using the y(A) constants estimated as in Figure 3. The dotted lines were computed by equation (12) that was calculated for 
10,000 sequential values of w in [0,1]. To this end, we replaced /i = 0 with /i= 10~ 6 in order to use equation (12) to plot the dashed red line. Initial 
conditions for the simulations: for the majority of the simulations we have simply used iv = 0.5 as initial condition for the plastic synaptic weight. In 
order to show the bi-stable solutions in panels (A2, B2, F1), for ;i = 0,0. 01,0.02 and a= 1.01,1.03, ... 1.19 we ran two simulations one with initial 
condition n =0 and another with initial condition w= 1. (A0-F0) are the STDP kernels (as in Figure 1) used in the simulations. (A1-F1) results for the 
inhibitory synapse simulations. (A2-F2) results for the excitatory synapse simulations. 
doi:10.1371/journal.pone.0101109.g005 
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Thus, the condition for instability of the ^-invariant solution is: 



dw 



(w*) = kr pre r post K{ 1 - w*)"(l + w*X + ) 



1 — w* 



l+w*X+ w* l+w*X_ 



(16) 



Yielding: 



sign 



d(w) 



dw 



(W) 



■■ sign 



dw* 
da 



,Vw*JT + >-l,vv*<l,a>0 (17) 



Hence, for w*X + > — 1, the fixed point solution in Figure 5 is 
stable in segments with negative slope, and unstable in segments 
with positive slope. Note that in our simulation setup \X + | < 1 (cf. 
Figure 4); thus, the condition w*X + > — 1 holds for all values of 8 
in our case. 

Revisiting the different scenarios depicted in Figure 5, we note 
the existence of two qualitatively different behaviors; namely, one 
that can only show mono-stability (A, C, and F) and the other has 
the potential for bi-stability (in panels B, D, and E). We use this 
behavior to classify the different STDP temporal kernels that are 
parameterized by the single variable 8. We shall term "class-I 
temporal kernels" the temporal kernels such that w* is mono- 
stable for all a>0; /xe[0,l]. We shall term "class-II temporal 
kernels" the temporal kernels such that w* is bi-stable for some 
flE[0,fl c ) and some a(ju). Note that this classification depends on 
the type of synapse (which via y(A) together with 8 determine 
X + i_). In addition, we note the existence of a special solution at 
w* = 1 1 2 that is invariant to \i, and enables us to obtain a simple 
condition for this classification. In class-I kernels the derivative 
dw* /da. at w* = 1/2 is always negative, whereas in class-II models 
there is a critical value of fi below which the derivative changes its 
sign. 

The "(i-invariant" solution and the critical \i 

As w* = l/2->(w* /{\-w*)y = Wn, the solution of the fixed 
point equation, equation (12), at w* = w= 1/2 is /i-invariant. For a 
given STDP temporal kernel (6), i.e. a given set of { X + ,X_ } (see 
Figure 4; and note that X + /_ are also determined by the pre-post 
correlation structure via y(A)), the solution of w= 1/2 is obtained 
with: 



, 1+X+/2 r „, n 



(18) 



Substituting the /i-invariant solution, equation (18), into 
equation (15), yields 



sign 



sign 



d<w>, 



(a,w) 



1 + X+/2 \+X_/2 



1 1 

1/2 + 1/2 



(19) 



VX+>-2 



1 ( x+ 



x_ 



4\\+X + /2 \+X_/2 



(20) 



Thus, for ji c < 0 the /.(-invariant solution, w, is stable for all 
values of ,ue[0,l] and the STDP rule is class-I for that synapse. On 
the other hand, if \i c > 0 the STDP rule is class-II. This 
classification depends solely on the values of {X + ,X-}. In our 
simulation setup \X + /_ | < 1 (see Figure 4), thus the classification of 
the parameter combinations is simply determined by the sign of 
(X + — X-); i.e. the manifold that is determined by the condition 
{X + =X-} separates the parameter space (that characterizes the 
STDP rule and the synapse) between class-I and class-II. 

Bimodal distribution near a 

Figure 6 depicts (using numerical simulations with set of class-II 
parameters) the bifurcation plots for the learning dynamics for 
inhibitory (A, B) and excitatory (C, D) synapses. For inhibitory 
synapses the anti-Hebbian (6 = — 1) plasticity rules were chosen, 
and for the excitatory synapses, the Hebbian (8= 1). The panels 
show the resultant distribution of the synaptic weight color-coded 
after 21x101 of 5 hours of simulations for 2 1 values of the 
bifurcating argument (either ji or a) along the abscissa. In order to 
calculate the synaptic weight distribution for the set of parameters 
without the bias of initial conditions, 101 simulations were 
performed with different initial weight values evenly spaced from 
0 to 1 . The rationale for running the simulations for 5 hours each 
was to make sure that the learning dynamics had reached a steady 
state regime and the synaptic weight fluctuated around it for the 
entire trailing 2.5 simulation hours. During these trailing 2.5 
simulation hours, the synaptic weights were recorded at a 1 Hz 
sample rate. For the estimation of the weight distribution, all the 
samples from the 101 simulations (differing only by their initial 
conditions) were used with 40 evenly spread bins between 0 and 1 . 

As expected from the analysis, there was a bifurcation along the 
H dimension (top panels), in which above fJ, c the distribution was 
uni-modal whereas below fi c the distribution was bi-modal. Along 
the a dimension (bottom panels) the distribution resembled the 
theoretical (dashed) curves of Figure 5 (without the unstable 
segment of dw* /da>0). 



Symmetry and phase transition along 0 

The high degree of similarity between the simulation results for 
inhibitory and excitatory synapses (Figure 5) stems from the fact 
that they obey the same mean-field equation (11), albeit with a 
different set of parameters. Thus, an excitatory synapse, w exc , with 
a specific choice of parameters {a,fi,X + ,X- } obeys the exact 
same mean-field equation as (1 — vv; n h), where H^h is an inhi- 
bitory synapse with the transformed set of parameters 
1 (1+X+) -X. -X+ 



,/<, 



and a somewhat different 



a (\ + X_)"~'\+X_ , \+X +- 
learning constant (note that X + j_ are positive for excitatory 
synapses and negative for inhibitory ones, see Figure 4). 

This symmetry is illustrated for different STDP temporal kernels 
in Figure 7, where the mean field fixed point, w* , is plotted as a 
function of a for different values of 8 (color coded) at ^~0. The 
different 8 were chosen around 8 C which is defined by the 
condition X + = X_ (see Figure 4) to display the phase transition 
from class-I to class-II along this parameter. Coincidentally, in our 
simulations and the chosen model (equations (7) and (8)), this 
specific 8 C was almost the same for excitatory and inhibitory 
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Figure 6. Bifurcation plots along the two parameters (,// z) of the weight dependent STDP component, / ( in (see equations (3) and 
(4)) near a of equation (18). Panels display the synaptic weight distribution (color coded) for the various parameter setups: (A) Inhibitory synapse 
with anti-Hebbian (0= — 1, see also Figured 1F) rule, with fixed a =1.07 and varied (B) Inhibitory synapse with anti-Hebbian (0= — 1, see also 
Figured 1 F) rule, with fixed \i = 0.025 and varied a. (C) Excitatory synapse with Hebbian (0=1, see also Figured 1B) rule, with fixed a = 1.13 and varied 
/(. (D) Excitatory synapse with Hebbian (0= 1, see also Figured 1B) rule, with fixed /i = 0.05 and varied a. The dashed white line marks fi c in A and B, 
and a in C and D. 

doi:10.1371/journal.pone.0101109.g006 

synapses; i.e. for both synapses 9 C > — 0.22 and 6 C < —0.21 (see 
Figure 4). Under these conditions, for an excitatory synapse, 
9< —0.22 defines the class-I kernels, and 6> —0.21 the class-II, 
whereas for an inhibitory synapse, 9< —0.22 defines the class-II 
kernels, and 8> —0.21 the class-I. 

Discussion 

The computational role of the temporal kernel of STDP has 
been studied in the past. Cateau and Fukai [8] provided a robust 
Fokker-Planck derivation and analyzed the effects of the structure 
of the STDP temporal kernel. However, their analysis focused on 
excitatory synapses and the additive learning rule (jl = 0). Previous 
studies have linked the Hebbian STDP of inhibition with negative 
feedback which acts as a homeostatic mechanism that balances the 
excitatory input to the postsynaptic cell [13,14]. Positive feedback 
and bi-stability of STDP dynamics have been reported only for 
excitation, and linked to sensitivity to the input correlation 
structure [6,7]. Here it was shown that the STDP of both 
excitation and inhibition can produce either positive or negative 
feedback depending on the parameters of the STDP model. Thus, 
for example, it was reported that both a temporally asymmetric 
Hebbian STDP (0=1) and a temporally symmetric learning rule 
(0 = 0) for inhibitory synapses generate negative feedback [13,14]. 
These reports are in-line with our finding that the critical 6 for 
transition from negative to positive feedback for inhibition is 
negative (6 c x —0.2). 

In general, STDP dynamics of single synapses was classified 
here into two distinct types. With class-I temporal kernels, the 




1.005 



Figure 7. Fixed point solution, u ', of the mean field approx- 
imation, (plotted using equation (12)), as a function of a, at 
,u«0 is shown for different values of 0 (color coded). Using pi>0 
yields continuity of the curves at the extreme values (w* = 1 and w* =0), 
which makes the picture clearer. On the other hand as the value of \i 
increases the unstable regime of a gets smaller and the resolution for 0 
steps plotted should decrease. Thus, to plot these lines, we used 
fi= 10~ 4 which is sufficiently close to 0 to illustrate the phase transition 
with high accuracy in 0. (A) Excitatory synapse. (B) Inhibitory synapse 
doi:10.l371/journal.pone.0101109.g007 
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dynamics is characterized by a negative feedback and has a single 
stable fixed point. In contrast, class-II temporal kernels are 
characterized by a sub-parameter regime in which the system is bi- 
stable (has positive feedback), and another sub-parameter regime 
with negative feedback. However, the mechanism that generates 
the negative feedback, (i.e., the stabilizing mechanism) in the two 
classes is different in nature. Whereas in class-I the negative 
feedback is governed by the convolution of the pre-post 
correlations with the temporal kernel, (i.e. the mean field constants 
X + /_, similar to the homeostatic mechanism in [13]), in class-II, 
the stabilizing mechanism is the non-linear weight dependent 
STDP component, /+ (w). Hence, there is no reason a-priori to 
assume that the negative feedback in class-II should act as a 
homeostatic mechanism. 

We found that there is no qualitative difference between the 
STDP of excitatory and inhibitory synapses and that both can 
exhibit class-I and class-II dynamics. Moreover, there is an exact 
symmetry between the excitatory and inhibitory STDP under a 
specific mapping of the parameters {cx,fi,X + ,X_ }. This symmetry 
results from the fact that the mean-field dynamics depend solely on 
the global mean field constants X + /_. It is important to note that 
although neural dynamics is rich and diverse, due to the separation 
of time scales in our problem, the STDP dynamics only depends 
on these fine details via the global mean field constants X + /_. 

Certain extensions to our work can be easily implemented into 
our model without altering the formalism. For example, empirical 
studies report different time constants for depression and 
potentiation, e.g. [1]. However, although in our simulations we 
used identical time constants at |0| = 1, for |f|<l the depression 
time constant is larger than the potentiation time constant in our 
simulations. Moreover, our analytical theory depends on the time 
constants only via X + /_. Consequently, changing time constants 
or any other manipulation to the temporal kernel can be 
incorporated into our mean-field theory by modifying X + 1 _ . 
Similarly, assuming separation of time-scales between short term 
and long term plasticity, the effect of short term plasticity can be 
incorporated by modifying X + /_ accordingly. 

STDP has also been reported to vary with the dendritic 
location, e.g. [18,25]. For a single synapse this effect can also be 
modeled by a modification of the parameters X + ; _ . However, the 
importance of the dendritic dependence of STDP may reside in 
the interaction with other plastic synapses along different locations 
on the dendrite. Network dynamics of a 'population' of plastic 
synapses is beyond the scope of the current paper and will be 
addressed elsewhere. 

In our model we assumed that the contribution of different 
"STDP events" (i.e., pre-post spike pairs) to the plastic synapse are 
summed linearly over all pairs of pre and post spikes, see e.g. 
equation (21). However, empirical findings suggest that this 
assumption is a mere simplification, and that STDP depends on 
pairing frequency as well as triplets of spike time and bursts of 
activity, e.g. [3,26-30]. The computational implications of these 
and other non-linear interaction of spike pairs in the learning rule, 
as well as the incorporation of non-trivial temporal structure into 
the correlations of the pre-synaptic inputs to the cell are beyond 
scope of the current paper. 

Empirical studies have reported a high variability of STDP 
temporal kernels over different brain regions, locations on the 
dendrite and experimental conditions, e.g., [1,12,15,17-19]. Here 
we represented the STDP rule as the sum of two separate 
processes, one for potentiation and one for depression with an 
additional parameter, 9, that allows us to continuously modify the 
temporal kernel and qualitatively obtain a wide spectrum of 



reported data. Representation of STDP by two processes has been 
suggested in the past. Graupner and Brunei [31], for example, 
proposed a model for synaptic plasticity in which the two processes 
(long term potentiation and depression) are controlled by calcium 
level. Thus, in their model the control parameter is a dynamical 
variable that may alter the plasticity rule in response to varying 
conditions. In our work, however, we did not model the dynamics 
of 9. Moreover, we assumed that 9 remains constant during 
timescales that are relevant for synaptic plasticity. It is, neverthe- 
less, tempting to speculate on a metaplasticity process [32,33] in 
which the temporal structure of the STDP rule is not hard wired 
and can be controlled and modified by the central nervous system. 
Thus, in addition to controlling the learning rate, X, or the relative 
strength of potentiation-depression, a, a metaplasticity rule may 
affect the learning process by modifying the degree of 'Hebbia- 
nitty', 9. Such a hypothesis, if true, may account for the wide range 
of STDP kernels reported in the experimental literature. How can 
such a hypothesis be probed? One option for addressing this issue 
is to try and characterize 9 during different time points and study 
its dynamics. One would expect to find that 9 (for excitatory 
synapses) decreases with time in cases where the neural network 
has been reported to becomes less sensitive to its input statistics, for 
example during developmental changes. 

Methods 

"Mean field" Fokker-Planck approach for the learning 
dynamics 

From the synaptic update rule, equation (5), changes in the 
synaptic weight, w, at time t, result from either pre or post synaptic 
firing at time t, affecting both the depression and potentiation 
branches (functions) of the adaptation rule. Thus: 



w(t + At)-w(t) 

-y.Mp(5,E[y+Ar]) 



dt'p pre (t')K + (t-t') 

-co 

/ 

f dt'p posl {t')K + {t'-t) 
L (21) 

dt'p pre (t')K_(t-t') 

— CO 
/ 

| dt>p posl {t')K_{t'-t) 



— CO 

where p pre / post (t) = E; A ( ? ~ <? re/part ) is the firing of the pre/ 

post synaptic cell, as represented by a train of delta functions at 
the neuron's spike times, with { tf } ^_ j being the spike times, 

and p(^' pj [^"' U e[t,t + At)^j is 1 if there was a pre/post synaptic 

spike respectively at the specified time interval [t,t + At) and 0 
otherwise. 

Taking the short times limit, At^O: hm w ' fe . = 

Ppre/posM, yields: 
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w(t)=Xf + (w)C + (t)-Xf_(w)C-(t) 



C + '-(t)= f p pre (f)p post (t)K +/ _(t-t')df 

J — CO 



+ 



i P P Jt)p po Jt')K +/ -(t'-t)dt' 

J — CO 



(22) 



(23) 



Assuming the learning process is performed on a much slower 
time scale than the neuronal dynamics [34], the STDP dynamics 
samples the pre-post correlations, p pre (t) P pos t(t) > over l° n S 
periods in which the synaptic weight, W, is relatively constant. 

Using this separation of time scales in the limit of X — >, we can 
approximate C + l~ (t) by their time average over period T« 1//L 
This is the mean-field Fokker-Planck approach to approximating 
the stochastic dynamics of W. Integration of equation (22) over 
time yields: 



lim 



L + l~= lim 



w(t)dt 
■ IT 



(24) 



dA 



dt 



+ 



f „ d A f T ^fP P re(t)p post (A+t)K +/ „(A)) (25) 



Assuming we can replace the time averaging of the pre-post 
correlation with its statistical average for sufficiently large T, 



r(A) = <p (0)/w(A)> = lim 



Ppre 

_ T 2T 



(26) 



we can substitute equation (26) into equations (25) and obtain the 
mean field Fokker-Planck equation for the process: 



<w>(w) = Xf + (w) m-A)K + ( — A) + T(A)K + (A)]dA 

J — CO 

,0 

-Xf-(w)\ [T(-A)K_(-A) + r(A)K_(A)}dA (27) 

J — GO 

/■CO PCO 

= ;/+(w) T{A)K + {A)dA-kf_{w) T(A)K_{A)dA 

J — 00 J — 00 



Details of the numerical simulations 

Online supporting information. This manuscript is ac- 
companied by a complete software package that was used 
throughout the study. This package is a Madab set of scripts 
and utilities that includes all the numerical simulations that were 
used to produce the figures in this manuscript. It also contains all 
the scripts that generated the figures. 

The leaky integrate-and-fire model. The learning dynam- 
ics of equation (5) was simulated by a single postsynaptic integrate- 



and-fire cell. As in our previous work [13] the dynamics of the 
membrane potential of the postsynaptic cell, V(t), obeys: 



dV _ {V rest -V) 

' dt ' 



+g E {E E -V)+g I (E I -V) (28) 



where C m =200pF is the membrane capacitance, R m = 100MQ is 
the membrane resistance, the resting potential is V rest = — 70mV, 
and the reversal potentials are Ee = OmV and Ej = — 70mV. An 
action potential is generated once the membrane potential crosses 
the firing threshold Vth = — SAm V, after which the membrane 
potential is reset to the resting potential without a refractory 
period. The synaptic conductances, gE and gj, are a superposition 
of all the synaptic contributions, i.e., each synaptic input is 
convolved with an oc-shaped kernel (that models the filtering 
nature of the synaptic response) amplified by its synaptic weight 
and then summed. The terms gE and gj are thus given by: 

^m=^e(^me H] + r ^/ tjr ) (29) 

where X stands for Excitation or Inhibition, Nx is the number of 
synapses, [t] + =max(t,0) is the dimensionless time value (in 

seconds), and | ? /|. are the spike times of synapse (. For the 

temporal characteristic of the a-shape response we chose to use 
T£ = T/ = 5 ms, and for the conductance coefficient g x our 
constant is scaled by Nx as elaborated below. 

In order to estimate the postsynaptic membrane potential in 
equation (28), the software performs the integration of the synaptic 
and leak currents using the Euler method with a At=\ms step 
size. The rationale for using such a low resolution step size and its 
verification are discussed below. 

Modeling presynaptic activity. Throughout the simula- 
tions in this work, presynaptic activities were modeled by an 
independent homogeneous Poisson processes, with stationary 
mean firing rate r pre = 10 spikes/ s. To this end, each of the inputs 
was approximated by a Bernoulli process generating binary 
vectors defined over discrete time bins of At= Ims. These vectors 
were then filtered using a discrete convolution oc-shaped kernel (as 
defined above) with a limited length of I0x x (after which this 
kernel function is zero for all practical purposes). In all simulations 
we used: N E = 120,iV/=40. 

Conductance constants. In order to be compatible with 
previous studies; e.g., [7,13], and to have simulations that are 
executed with a robust and generic software package accompa- 
nying this manuscript as File SI, we scaled the synaptic 
conductance inversely to the number of synaptic inputs in our 
simulations. We used the following scaling formula g x = g x Sx, 
with: g§ = 30nS, S E = W00/N E , gf = 50nS and S/ = 400/7Y/, 
where Ne,Ni are the number of excitatory and inhibitory 
presynaptic inputs, respectively. 

The learning rate. The simulations of the STDP process 
were carried out to obtain the asymptotic weight distribution of the 
plastic synapse. Convergence to the asymptotic region was 
accelerated by manipulating the learning rate constant k of 
equation (1). The software code was designed to support a given 
vector of X for each minute of the simulation. Specifically 
we used the following formula to generate this vector: 

10~ 3 \\ +9(1 -?) 10 j , where fe(0,l], is the ratio between the 

minute iteration time and the entire simulation time. Examining 
the behavior of this function shows that it starts from a value of 
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10 2 and decays significantly fast, leaving the trailing 70% of the 
simulation time with more or less the same learning rate of about 

io- 3 . 

Postsynaptic spike time accuracy vs. simulation step size 
resolution. Figure 5 shows the remarkable match between the 
fixed point solution (w*) of equation (12), and the asymptotic 
synaptic weight (wq) of the simulations; the regression coefficient 
on the entire set {wo,w*} in all the panels is 1 + 5 - 1 0 4 with 
R 2 > 0.999 when using an integration step of size Ims. Tests of this 
kind were performed on simulations using integration steps 
ranging from 0.1ms to Ims in two calculation modes (see below), 
and it was found that higher resolution provides a better match to 
the analytical solution. However, the key feature that contributes 
to this high degree of similarity between the analysis and the 
simulations (more than an order of magnitude for the error term 
1 — R 2 ) was the definition of the spike times of the postsynaptic cell 
rather than a 10 x decrease of the integration step size. 

The spike times of an integrate and fire neuron are defined as 
the times in which its membrane potential crossed the firing 
threshold, t* . However, in the numerical simulations we used 
discrete times, {wAfbinl^Lo- 1° previous work we define the time 
of the post-synaptic firing by the last discrete time preced- 
ing the threshold-crossing time to: t post = nAt such that 
nAt< f < (n+ l)At. This choice may change the causal order of 
pre-post firing (from pre before post to simultaneous firing) at time 
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